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TThis  report  describes  a set  of  computer  programs  which  calculates  the 
inviscid  two-dimensional  expansion  of  gaseous  exhaust  mixtures  through  rocket 
nozzles.  Equilibrium  chemistry  in  the  combustion  chamber  and  as  far  as  the 
nozzle  throat  is  assumed,  whilst  finite  rate  or  frozen  chemistry  may  be  applied 
from  the  throat  onwards.  The  transonic  flow  field  is  calculated  by  the  method 
of  streamline  curvature  and  the  flow  in  the  supersonic  region  of  the  nozzle  is 
solved  by  the  method  of  characteristics. 


The  theory  on  which  the  programs  are  based  is  summarised  and  a comprehen- 
sive user's  guide  is  included.  A sample  calculation  is  presented. 
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The  computer  programs  described  provide  means  of  calculating  the  steady, 
inviscid  two-dimensional  expansion  of  a gaseous  mixture  through  an  axisymmetric 
rocket  exhaust  nozzle  (Fig.  1),  allowing  for  the  effects  of  chemical  reaction 
rates  during  expansion  of  the  exhaust  gas. 

The  set  of  programs  consists  of: 

(a)  ODNOZ  (one-dimensional  nozzle  expansion),  which  permits  the  rapid  calcula- 
tion of  the  equilibrium  temperature  and  composition  inside  a rocket  motor 
combustion  chamber,  and  provides  an  axial  profile  of  gas  composition, 
temperature  and  pressure  along  the  nozzle  with  the  assumption  of  equilibrium 
chemistry  (infinite  reaction  rates).  Tables  of  flow  properties  (x,T,p,m,M) 
and  species  mole  fractions  X.  are  constructed. 

(b)  TRANOZ  (transonic  nozzle  flow),  which  calculates  an  average  expansion 
coefficient  y from  the  above  tables,  and  assuming  constant  stagnation 
conditions  calculates  the  axisymmetric,  transonic  flow  field  by  the  method 
of  streamline  curvature.  An  initial  constant  property  data  line  is 
constructed  for  input  to 

(c)  SUPERNOZ  (supersonic  nozzle  program),  which  calculates  the  flow  in  the 

supersonic  region  of  the  nozzle  by  the  method  of  characteristics,  with 

either  a frozen  chemical  composition  (zero  reaction  rates)  or  taking  into 

account  the  finite,  non-zero  rates  of  chemical  reactions.  It  provides 

accurate  exit  plane  conditions  which  may  be  used  as  input  to  an  exhaust 

plume  analysis  program,  such  as  the  Rocket  Exhaust  Plume  program  (REP) ^ or 

2 

the  Base  Flow  program  (BAFL)  . 

The  sequence  of  calculations  is  illustrated  in  Fig.  1. 


i 


The  programs  already  represent  a significant  advance  over  those  currently 

used  at  PERME  for  calculating  the  flow  from  the  combustion  chamber  to  the  nozzle 

exit  plane^,  in  that  they  can  treat  two-dimensional  flow,  and  they  are  to  be 

further  enhanced  by  including  the  non-equilibrium  effects  of  gas/particle  flow. 

These  have  already  been  incorporated  into  the  transonic  flow  program  wi thout 
I* 

chemistry  , and  the  full  two-dimensional,  two-phase  flow  treatment  should 
shortly  be  available. 

This  report  serves  primarily  as  a program  user's  guide,  but  also  contains  a 
summary  of  the  theoretical  basis  of  the  programs.  The  governing  equations  are 
included  and  the  method  of  solution  for  each  of  the  three  programs  is  described. 
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The  data  and  results  of  a sample  calculation  are  presented.  The  results  are 
all  given  in  SI  units,  and  the  working  units  for  the  fluid  dynamic  equations 
are  SI.  Ir  the  chemical  equations  most  of  the  variables  are  non-dimensional. 
The  programs  are  written  in  FORTRAN  IV,  permitting  their  use  on  all  standard 
computers  and  are  designed  for  easy  use,  requiring  only  specification  of  the 
propellant  composition,  pressure  and  enthalpy,  nozzle  goemetry,  thermodynamic 
data  and  the  chemical  reaction  scheme  with  rate  coefficients. 


2 ONE-DIMENSIONAL  EXPANSION  WITH  EQUILIBRIUM  CHEMISTRY 

The  program  ODNOZ  is  a modified  version  of  NEWFEC^,  with  the  addition  of  a 
cubic  spline  fitting  routine,  which  permits  a contoured  nozzle  from  the  throat 
onwards.  Most  of  the  theory  and  methods  used  in  the  program  are  taken  from 
Ref.  5,  any  differences  in  the  chemical  treatment  being  described  in  Ref.  3- 


2.1  Chamber  calculation 


The  gaseous  propellant  mixture  is  made  up  of  s chemical  species,  of  which 
c are  either  chemical  elements  (H,0,N,  etc.),  or  linearly  independent  combina- 
tions of  the  chemical  elements  (Hj^.Nj,  etc.),  and  s-c  are  dependent  species 
formed  by  simple  combinations  of  the  c independent  species.  The  system  is 
specified  through  the  matrix  A , in  which  the  element  represents  the 

number  of  moles  of  independent  species  j per  mole  of  dependent  species  i 
The  number  of  moles  of  the  fundamental  species  per  unit  mass  of  the  unreacted 
mixture  q(  , j = 1 , 2,  ...  c , must  be  given.  The  conservation  of  chemical 
elements  and  global  mass  conservation  are  expressed  by  the  equations 


X.  = q. 


i-Z 

i=c+1  L-  k=1 


- I 


X.  j-1,  2,  ...  c (2.1) 


where  X.  , i = 1 , 2,  ...  s is  the  mole  fraction  of  species  i and 


q;  = 


q • 


J 


j=1 


At  chemical  equilibrium  the  dependent  species  concentrations  are  related  to  the 
concen t ra t ions  of  the  components  through  the  equilibrium  constants  in  the 
equation: 
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i“C+1 , . . . s 


(2.2) 


where  is  the  mixture  pressure  in  atmospheres.  The  factor 


exp 


>1 


o 

o 

p . 

V, 

a — L 

i 

ij  T 

T 

J n 

n 

is  the  equilibrium  constant  based  on  partial  pressure  for  species  i 


n 


i s 


the  mixture  temperature  T , normalised  by  the  combustion  chamber  temperature 

o 


T , so  that  in  the  chamber  T = 1 , and  p.  is  the  normalised  chemical 

o n i 


potential  of  species  i 


Substitution  from  equation  (2.2)  into  equation  (2.1)  for  the  X. 


i=c+1,  ...  s , reduces  equation  (2.1)  to  a set  of  c simultaneous  equations 
which  may  be  solved  iteratively  for  the  equilibrium  concentrations  of  the 
independent  species  at  a qiven  pressure  P and  temperature  T . The  chamber 

3 

temperature  T is  not  known  initially  and  so  the  specific  enthalpy  H 


o ' ■ ■ o 

which  may  be  readily  determined  from  the  heats  of  formation  of  the  constituent 


compounds,  is  specified  and  the  solution  is  found  by  iteratively  varying  the 
temperature  until  the  mixture  enthalpy  H matches  the  specific  enthalpy. 


2 . 2 Equilibrium  expansion 

The  equations  of  continuity  and  conservation  of  momentum  and  energy  for  an 
inviscid  gas  are 


pu  A = constant 


u du  + — dP 

P 


(2.3) 

(2.4) 


H + |u2  « H 

2 o 


constant 


(2.5) 


where  u is  the  velocity,  p is  the  density  and  A is  the  cross-sectional 
area  of  the  nozzle  normalised  by  the  throat  area. 


T ^ 
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The  equation  of  state  for  an  ideal  gas  is 


p T R 


where  m is  the  molecular  weight  and 


is  the  universal  gas  constant. 


The  energy  and  momentum  equations  may  be  used  to  show  that  the  inviscid 
expansion  of  a flow  in  chemical  equilibrium  is  an  isentropic  process.  The 
molar  entropy  of  a chemically  reacting  gas  mixture  is  given  by 


s s 

S - - I > X.  s°.  - > X.  In  X.  - In  P 

m / Z_  J J L—  J J a\ 


(2.7) 


where  s.  is  the  molar  entropy  of  species  j at  1 atmosphere.  Equations  (2.1) 
and  (2.7)  now  give  a set  of  c+1  simultaneous  equations  for  determining  the 
composition  and  pressure  of  the  mixture  at  a given  temperature.  The  equilibrium 
solution  is  obtained  at  successive  temperature  intervals,  computing  the  flux, 
pu  , at  each  point  from  equations  (2.5)  and  (2.6).  The  maximum  flux  and  the 
corresponding  temperature  determine  the  throat  conditions  and  the  critical  mass 
flow  for  the  equilibrium  expansion.  After  the  mass  flow  has  thus  been  found 
the  chemical  composition  and  aerodynamic  variables  at  given  axial  stai  ons  may 
be  calculated  from  equations  (2.1)  to  (2.7),  and  tables  of  flow  properties 
(x,T,p,m,M)  and  species  role  fractions  X.  , i-1,2,  ...  s are  constructed  for 
input  to  the  next  stage  of  the  calculation.  The  equilibrium  exit  plane  condi- 
tions are  also  given. 

3 TWO-DIMENSIONAL  TRANSONIC  FLOW 

In  the  solution  of  the  two-dimensional  equations  for  flow  through  a 
convergent/divergent  nozzle  with  a transonic  region,  many  mathematical  difficul- 
ties are  encountered  in  the  region  of  transition  from  subsonic  to  supersonic 
flow.  Although  the  transition  is  smooth  the  properties  in  the  two  regions  are 
quite  different.  In  the  subsonic  region  the  governing  equations  are  elliptic, 
whereas  in  the  supersonic  region  they  are  hyperbolic.  The  supersonic  flow  field 
may  be  solved  by  the  method  of  characteristics,  but  accurate  initial  line  data 
are  needed  to  start  the  calculation.  Since  the  flow  near  the  nozzle  throat  may 
be  far  from  one-dimensional,  the  sonic  velocity  at  the  wall  being  reached  well 
upstream  of  the  throat,  while  that  at  the  axis  is  reached  downstream  of  the 
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thro.it  (see,  for  example.  Fig.  2),  a method  of  solution  which  describes  the 
properties  of  both  subsonic  and  supersonic  flow  segments  must  be  devised.  A 
further  difficulty  is  that  the  mass  flux  is  not  known  and  must  be  calculated. 
When  chemical  reactions  are  taken  into  account  the  problem  becomes  even  more 
complex  and  it  is  thus  necessary  to  construct  by  approximate  means  an  initial 
line  suitable  for  starting  the  method  of  characters  tics  solution  of  the  flow 

* field  in  the  supersonic  domain. 

v In  TRANOZ,  an  effective  expansion  coefficient  is  computed  by 

i 

Y “ 

where  L refers  to  the  last  entry  in  the  table  constructed  by  ODNOZ  and  o 
denotes  the  chamber  value.  A mean  molecular  weight  is  found  f.om 

i - „ , 

Z o L 

i 

• and  the  stagnation  conditions  are  assumed  to  remain  constant  at  the  chamber 

val ues . 

I 

] » 

There  are  several  ways  of  solving  the  equations  of  transonic,  axisymmetric 
flow,  but  most  are  restricted  to  nozzles  of  a particular  shape  and  with  a throat 
wall  radius  of  curvature  greater  than  (at  least)  the  throat  radius.  The  method 
adopted  in  TRANOZ  is  based  on  that  described  by  Stow^,  which  was  developed  from 
the  "streamline  curvature"  method  used  for  subsonic  flows  in  turbomachines  (see, 
for  example.  Marsh  ).  It  is  an  iterative  procedure  which  may  be  applied  to 
nozzles  of  general  geometry  and  which  is  convergent  and  stable,  even  for  nozzles 
with  a throat  wall  radius  of  curvature  smaller  than  the  throat  radius. 

When  the  transonic  flow  field  has  been  calculated,  the  wall  Mach  number  at 
a specified  point  is  found  and  a table  of  (x,r,0)  values  corresponding  to  this 
value  of  M on  each  streamline  is  constructed,  to  start  the  supersonic  solution 
(see  Fig.  2).  The  point  chosen  is  usually  the  throat  wall  point,  since  this 
permits  a finite  rate  chemical  treatment  as  far  upstream  as  possible.  The 
initial  line  thus  constructed  is  an  approximation  to  the  constant  pressure 
surface  for  a Mach  number  of  slightly  greater  than  unity,  along  which  it  is 
reasonable  to  assume  that  a constant  expansion  coefficient  can  be  used.  The 
other  gas  properties  corresponding  to  this  Mach  number  are  obtained  by  inter- 
polation from  the  tables  constructed  by  ODNOZ. 


In  (P,/P  ) 
L o 


In  Ip.  /pi 

L O 
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3. I Method  of  streamline  curvature 

The  equations  governing  the  steady  ax  i symme  t r i c motion  of  an  inviscid 
compressible  fluid  are,  in  the  co-ordinates  of  Fig.  3,  where  x and  r are 
normalised  by  the  throat  radius  r*  : 


cont inui ty 


h (rpu)  + fr  (rpv) 


conservation  of  nomen  turn 


(3.1) 


i • 


3u  3u 

u — ♦ V — 
3x  3r 


3v  3v 

u — ♦ v — 
3x  3r 


_ J_  3P 

P 3x 

1 if. 

p 3r 


Equation  (3.1)  implies  the  existence  of  a stream  function  4>  satisfying 


34/ 

— - rPu 


- rp  v 


or 


(x,  r) 


J (rpu) 


x=cons  tant 


dr  + constant 


1/2 


(3.2) 

(3.3) 


(3.4) 


By  introducing  the  speed  of  sound,  defined  by  c = (3P/3p)  , equations  (3-1) 

to  13.3)  may  be  combined  to  obtain  the  radial  equilibrium  equation,  relating  the 
static  pressure  field  to  the  streamline  geometry.  Detailed  derivations  are 
given  in  Refs  ^ and  6,  but  only  the  final  equation  is  reproduced  here*: 


f 1 u2  2 1 

1 ( 3P 

1 

2 

! - u2  1 

[ 3 tan  0 

1 - — 55-  sec  e 

*■  c 

7 i 3r 

lx  " 

- u 

17' 

l 1 

2.  , f 3 tan  0 

♦ u tanfl  , 

l 3r 


2 , 5 

u tan 


(3.5) 


/"/Subscripts  x , 4>  imply  differentiation  along  a line  of  constant  x or  4' 
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The  term  — — — is  referred  to  as  the  "curvature"  of  a streamline. 

i 3x  )* 

3.2  Solution  procedure 

Equation  (3.5)  is  solved  iteratively  along  axial  stations  in  the  flow. 

Each  axial  station  is  divided  into  a number  of  radial  mesh  points  (see  Fig.  3) 
and  the  equation  is  applied  at  points  mid-way  between  the  mesh  points.  It  is 
assumed  that  an  estimate  for  the  flow  pattern  is  available  from  a previous 
■ iteration,  so  that  equation  (3.5)  becomes  an  ordinary  differential  equation  in 

the  static  pressure. 

To  begin  the  iterative  procedure  an  approxi mat  ion  is  made  to  the  mid-radial 
pressure  and  the  shape  of  the  streamlines  f rom  one-dimens ional  considerations, 
and  equation  (3.5)  is  integrated  radially  to  obtain  the  pressure  profile.  The 
mass  flux  is  determined  from  equation  (3.M  and  the  mid-radial  pressure  is 
varied  until  the  calculated  mass  flux  agrees  with  the  true  mass  flux.  All  the 
axial  stations  are  solved  in  this  way,  and  the  other  flow  variables  are  evaluated 
from  isentropic  relationships.  A new  approximation  to  the  streamline  geometry 
is  made  and  the  calculation  is  repeated  until  the  solution  converges. 

In  the  case  of  choked  internal  flow  the  mass  fl>tx  must  be  determined. 

This  is  calculated  at  the  throat  in  each  iteration  by  choosing  the  mid-radial 
pressure  which  maximises  the  mass  flow  through  the  throat.  This  is  equivalent, 
in  one-dimensional  flow,  to  finding  the  sonic  point.  The  calculation  is  then 
progressed  upstream  and  downstream  of  the  throat  using  this  mass  flux. 

At  axial  stations  upstream  and  downstream  of  the  throat  two  possible 
isentropic  solutions  exist.  These  correspond  to  the  subsonic  and  supersonic 
solutions  of  one-dimensional  flow,  although  in  the  two-dimensional  case  the  flow 
can  be  mixed  in  the  radial  direction.  Upstream  of  the  throat  the  equivalent  of 
the  one-dimensional  subsonic  solution  is  chosen  and  for  rocket  nozzles,  where 
the  back  pressure  is  such  that  the  flow  becomes  supersonic  downstream  of  the 
throat,  the  equivalent  of  the  one-dimensional  supersonic  solution  is  chosen  at 
stations  downstream  of  the  throat. 

I In  practice  the  iterative  procedure  is  conditionally  stable  and  relaxed 

values  of  the  calculated  slope  and  "curvature"  must  be  used  to  ensure  convergence. 
I • For  example 
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(tant>)n  , •*  (tanu)n  ' ♦ relaxation  factor  x 

used  used 

1 (tanO)"  . - (tanO)n  ' | 

[ found  used  J 

where  (tan*})n  represents  the  nth  approximation  to  tan  6 and  the  subscripts 
"found"  and  "used"  refer  to  the  value  calculated  and  that  used  in  the  next 
iteration.  The  relaxation  factor  depends  upon  the  nozzle  geometry,  as  does  the 
number  of  steps  required,  convergence  being  assumed  when  the  change  in  the  flow 
parameters  between  successive  iterations  is  smaller  than  some  desired  degree  of 
accuracy . 

j 

k SUPERSONIC  METHOD  OF  CHARACTERISTICS  SOLUTION  WITH  NON-EQUILIBRIUM  CHEMISTRY 

The  computer  program  SUPERNOZ  calculates  the  inviscid  two-dimensional  frozen 
or  non-equi 1 ibrium  expansion  of  a gaseous  mixture  through  the  supersonic  region 
of  a rocket  exhaust  nozzle  by  the  method  of  characteristics.  Characteristics  are 
lines,  like  streamlines,  along  which  the  partial  derivatives  may  be  grouped 
together  in  such  a manner  that  they  form  ordinary  derivatives  in  the  direction 
of  the  charac te ri s t i cs . For  supersonic  flow  the  charac te r i s t 1 cs  are  real  and 
coincident  with  the  Mach  lines;  thus,  for  any  point  P in  the  flow  the  angle 
between  the  tangent  to  a characteristic  and  the  tangent  to  a streamline  is  the 
Mach  angle  n . For  every  such  point  P in  the  flow  there  are  two  character- 
istic lines,  referred  to  as  the  left  and  right  running  characteristics  (Fig.  k) . 

With  the  initial  data  line  constructed  by  TRANOZ  as  a starting  point,  the 
characteristic  equations  governing  the  fluid  dynamic  variables  are  integrated 

g 

via  an  explicit  integration  technique  , while  the  chemical  relaxation  equations 
are  integrated  using  an  implicit  integration  method^  to  ensure  numerical  stabil- 
ity in  near-equilibrium  flows.  In  explicit  numerical  solutions  the  solution  is 
calculated  step  by  step  from  the  given  differential  equation  and  known  initial 
conditions,  whereas  in  implicit  solutions,  the  unknown  values  are  bound  together 
by  a system  of  simultaneous  equations.  The  fluid  dynamic  equations  are  inte- 
grated explicitly  since  the  required  integration  step  sizes  are  not  stability 
limited,  and  an  implicit  formulation  of  the  equations  would  unnecessarily 
complicate  the  calculations. 

The  theory  on  which  the  solution  of  the  fluid  dynamic  equations  is  based  is 
taken  from  Ref.  8 and  the  integration  scheme  for  the  non-equilibrium  chemical 
equations  is  taken  from  the  one-dimensional  program  FiRAC^. 
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UNI  IMIUD 


4 • ' Ch. i i.u  tot  i s t it  i'i|u.i  t ions 

The  conservation  equations  governing  the  steady  axi symmel ri c Invlscid  flow 

g 

of  a react in«|  gas  mixture  ore.  In  characteristic  form  , 


dr 

— — - ton  0 

dx 


. V dP 

d — + — 

2 p 


dP 

yP 


dp 

t> 


cos  0 


dx 


) - I dP  dT 

Y P " T 


B 


cosO 


dx 


w . r* 

dc  . « — f; — dx 

j pvcosO 


(4.1) 


(4.2) 


(4.3) 


(4.4) 


(4.5) 


along  streamlines,  where  Cj  Is  the  species  mass  fraction,  Is  the  net 

species  production  rate  and  r*  is  the  throat  radius, 


4—  " cot  (0+u) 

dr 


1 

- c. 

’ f A - ilQS.  i 

F dr  - do 

i r | 

dP 

P 


along  left  running  characteristics,  and 
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(4.7) 


—r~  - ton  (o-a) 

dx 
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r 

dP 

P 


olon<|  right  running  characteristics,  where 
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M 
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Cp  - R/nt 


t 

F - com*  - slnO  cot  (0+a) 


( '» . M) 

(A. 12) 

(k.  13) 

(k.lk) 

(k.  15) 

(k.  16) 

(*».  17) 


a - — 1 

sina  cosn 


(‘•.18) 


H - cosO  tan  (0-o)  - slnO 


(k.  19) 


I ho  above  relationships  remain  determinant  provided  that  the  I low  is  always 
supersonic,  that  is  M > I 

I* . 2 Non-equ  i I i ht  i uni  t hem i s t vy. 

Ihe  non-equ 1 1 I hi  I urn  conditions  are  expressed  as  differential  forms  of  the 

, 5 

conservation  equations  : 
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0 k-1 ,2,  . ..  c 


(^ . 20) 
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(A. 21) 
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where  Is  the  concentration  variable  in  moles  per  gram,  defined  by 


X.  c. 


j m nu 


r is  the  number  of  reactions,  and 


0.  . - v!  . - v.  . 

■J  'J  ij 


r*  (.001  p) 
(100  u) 


s v,  , 


kf  TT  Yj  J moleg 
i j*>  1 J 


, , v Pi  s B.  . 
jp/£o)_TT  .j 

Ki  j-1  J 


where  v.  . , vl , are  the  stoichiometric  coefficients  of  species  j in  the 
ij  ij 

reaction  i on  the  reactants  and  products  sides  respectively,  is  the 

i 

forward  rate  coefficient  of  the  i th  reaction  and  K.  is  the  equilibrium 
constant  of  reaction  i . In  addition, 


V, 


•■-I*.;  ■ 

j-’ 

4 . 3 Numerical  integration  technique 

The  reacting  gas  characteristic  relationships  (equations  (A.l)  to 
are  usually  integrated  via  second  order  explicit  methods.  However,  for  near- 
equilibrium  chemistry,  the  explicit  form  of  equation  (4.5)  requires  impractic- 
ally  small  integration  step  sizes  to  ensure  stability.  Thus  the  fluid  dynamic 
equations  ((4.1)  to  (4.4)  and  (4.6)  to  (4.9)),  which  are  not  stability  limited 
to  small  step  lengths,  are  integrated  by  explicit  methods,  since  an  implicit 
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formulation  of  the  equations  vould  unnecessarily  complicate  the  calculations, 
and  the  chemical  concentrations  are  found  via  an  Implicit  scheme. 

4.3.1  Integration  ot  fluid  dynamic  equations 

Consider  the  flow  field  shown  in  Fig.  4,  where  point  4 is  a known  back  data 

point,  points  I and  2 are  known  points  located  on  the  right  and  left  running 

characterist ics  respectively,  and  point  3 is  the  point  to  be  calculated. 

Between  points  3 and  4 the  streamline  characteristic  relationships  are 

. 8 

integrated  as 


r^  + tan  j (0^  + 0^)  (x?  - x^) 


(4.22) 


\ w 2 2 

~|  V4  + 


Y-1 


1 - 


lzJL 

Y 


1/2 


(4.23) 


where 


In  (W 
In  (pj/p^) 


p4 


1 1 


3 ) 2 ^ y4  y3 


exp 


3 


cos0,.  cost). 


(*3  - x^) 


(4.24) 


T,  ^ 


V’  '3'' 


exp 


cosO^  cost)^ 


(x3  - X||) 


(4.25) 


Between  points  1 and  3 the  right  running  characteristic  relationships  are 
integrated  as 


r^  + tan 


(6j  + 0j  - - «3)  (x^  - Xj) 


(4.26) 
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« P(  exp 


' [l  (AlclHl  * ajg3h3>  <*3  ' *i> 


G.H , sinO,  G.H  j sinO, 

-!— ! — * -i-2 1 u,  - x.) 

r,  r3  3 I 


j (G,  + Gj)  (03  - 6,) 


If  point  3 is  an  axis  point,  and  0^  are  zero  and  the  indeterminate 

quantity  sinO^/r^  is  approxi mated  by 


r j + (x^  - x ^ ) tan61 


Between  points  2 and  3 the  left  running  characteristic  relationships  are 
integrated  as 


x^  * x2  + cot  T ^°2  + ^3  + a2  + a3^  ^r3  ” r2^ 


P^  = P0  exp  y + A3G3F3)  (r^  - r^) 

. r G,F_  slnO_  G,F-  sinO, 

. 1 2 2 2 3 3 5 / r - r ) 

7 l ^ + r^  r3  V 


j (g2  + g3)  (e3  - e2) 


If  point  2 is  an  axis  point  the  indeterminate  quantity  slnO^/^  is  that 
quantity  previously  calculated  for  the  axis  point  by  equation  (4.28). 

Equations  (4.27)  and  (4.301  can  be  combined  to  give 
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♦ y ( A ^ G 1 H ^ + A^G^H^)  ( x ^ - Xj) 


1 

7 


G)H1  sine 


j G3H3  s‘ne3  T 
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(*3  ‘ *j) 


(i*.  3D 


Since  the  above  inteqration  equations  involve  the  flow  properties  at  the 
unknown  points  3,  they  must  be  solved  by  iteration.  The  iterative  technique  Is 
described  in  Section  **.4. 

4.3.2  Inteqration  of  chemical  equations 

The  set  of  equations  qivinq  the  increments  Sy.  of  the  concentration 

3 1 

variables  y.  is: 
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(4.32) 


dy . 


where  f. 


i dx 


(5y | ) n+ 1 ” change  in  yj  in  a distance  (Ax)n+1  along  a streamline 


( A x ) . - x , “ x 

n+1  n+1  n 
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The  partial  derivatives  - — are  determined  numerically  from 


'^i  (Ypfjf  Y j + Y.)  “ f[  (Y»tT?f'«»  Y;  ” *;*•••»  Ye) 
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where  A.  is  some  small  fraction  of  v.  and  the  ordinary  derivatives  -r — 

J J dx 

are  found  by  solving  equations  (4.20)  and  ( 4 . 2 1 ) . Equation  (4.21)  cannot  be 
used  at  the  initial  data  line,  where  the  flow  is  in  equilibrium,  so  an  additional 
equation  is  obtained  for  starting  the  solution  by  differentiating  equation  (2.2), 
to  give 
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a,.  - 1 ] d-ln-i£l ,P°i  » 0 i-c+1,  ...  s 

i j J dx 


dT  din  ( t I ) 

where  the  derivatives  -r—  and  n , ' 0 1 ■ ° are  provided  along  the  initial 

dx  dx 

data  line  by  TRANOZ. 
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4 . 4 Solution  procedure 

The  characteristic  mesh  for  the  supersonic  solution  is  shown  in  Fig.  4. 

The  solid  dots  are  known  points  on  the  back  data  line,  and  the  squares  represent 
the  points  to  be  computed.  Point  2 is  the  point  on  the  back  data  line  next 
below  point  4,  and  the  unknown  point  3 lies  at  the  intersection  between  the 
streamline  through  point  4 and  the  left  running  characteristic  through  point  2. 
Point  1 is  at  the  intersection  between  the  right  running  characteristic  through 
point  3 and  the  left  running  charac ter i s t i c through  point  4.  Values  at  point  1 
are  determined  by  interpolation. 

The  iterative  solution  is  started  by  assuming  that  the  previously 
uncalculated  flow  variables  on  the  right  hand  sides  of  equations  (4. 22)  to  (4.31) 
are  equal  to  the  values  at  points  4,  2 or  1 depending  on  whether  the  integration 
is  along  a streamline  or  a left  or  right  running  characteristic,  respectively. 
Initial  estimates  of  the  locations  of  point  3 and,  except  for  wall  points, 
point  1 are  made  and  the  fluid  dynamic  variables  at  point  3 are  found  from 
equations  (4.22)  to  (4.31).  The  species  concentrations  at  point  3 are  calculated 
by  integration  of  the  chemical  equations,  and  the  other  properties  at  point  3 
are  found  from  equations  (4.10)  to  (4.19).  New  estimates  of  the  locations  of 
points  3 and  1 are  made,  and  the  solution  procedure  is  repeated,  always  using 
the  latest  computed  values  of  the  flow  properties  at  point  3.  until  r^,x^,0^, 
P3,P3’~3  3nC*  V3  conver9e  to  the  required  accuracy. 

The  first  point  to  be  computed  is  the  nozzle  wall  point.  Only  two  known 
data  points  are  required  here  since,  once  the  axial  location  of  point  3 has  been 
found,  r^  and  0^  may  be  determined  from  the  nozzle  geometry.  The  internal 
points  are  then  computed,  starting  with  the  point  adjacent  to  the  wall  and 
continuing  with  the  one  immediately  below  it,  until  the  nozzle  axis  is  reached. 
Again,  only  two  back  points  are  needed  here,  since  r^  and  0^  are  both  zero 
on  the  axis.  When  all  the  forward  points  have  been  calculated,  the  new  data 
line  becomes  the  back  data  line,  and  the  computation  is  progressed  downstream 
until  the  exit  plane  is  reached. 

4.5  Thrust  calculation 

When  the  exit  plane  conditions  have  been  found,  the  thrust  may  be  calculated. 
Consider  the  forces  on  a stationary  motor.  The  surrounding  atmosphere  exerts  a 
pressure  on  all  external  surfaces  and  at  the  nozzle  exit  plane  the  pressure 

is  that  of  the  exhaust  gas.  There  is  thus  a nett  forward  force  on  the  motor, 
which  is  equal  to  the  amount  by  which  the  gas  pressure  at  the  exit  exceeds  the 
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atmospheric  pressure,  integrated  over  the  nozzle  exit  plane.  The  rearward 
momentum  of  the  gas  stream  emerging  from  the  nozzle  produces  a reacting  forward 
force  proportional  to  the  momentum  flow  rate  of  the  gas.  Thus  the  total  forward 
force  on  the  rotor,  the  thrust  T , is  given  by 


L 


FIN 


j [rpu  + r (P  - PA)] 

1 o F 


d r 


IN 


5 USE  OF  THE  PROGRAMS 

The  three  computer  programs  are  interdependent  and,  with  the  exception  of 
ODNOZ,  which  may  be  used  separately  to  generate  one-dimensional  equilibrium  exit 
plane  conditions,  must  be  used  in  the  sequence  shown  in  Fig.  1.  They  are 
relatively  simple  to  use,  requiring  only  specification  of  the  propellant  composi- 
tion, pressure  and  enthalpy,  nozzle  geometry,  thermodynamic  data  and  the  chemical 
reaction  scheme  with  rate  coefficients.  The  thermodynani c data  and  reaction 
rate  coefficients  are  input  via  the  standard  data  files  THERMDATA  and  REACTIONS, 
which  are  searched  by  the  programs  to  find  the  thermodynamic  data  for  the  given 
species  and  the  rate  coefficients  for  the  given  reactions. 

The  nozzle  geometry  may  be  one  of  two  types,  featuring  either  a conical 
exit  section  (Fig.  5a)  or  a contoured  exit  section  (Fig.  5b).  If  the  nozzle  is 
contoured  a least-squares  approximation  is  made  to  a set  of  given  data  points  on 
the  nozzle  wall  (x.  , r.)  , with  weights  w.  , i-1,2,  ...,  m , by  the  method 

of  cubic  splines.  A cubic  spline  consists  of  a set  of  cubic  polynomial  arcs 
joined  end  to  end  with  continuity  in  value  and  in  first  and  second  derivatives 
at  each  join  (or  knot).  The  positions  of  the  interior  knots  are  chosen  by  the 
user.  The  algorithm  for  fitting  cubic  splines  to  a set  of  data  points  is  taken 
f rom  Ref.  9 . 

If  the  values  of  the  data  points  are  to  be  read  from  graphs,  sufficient 
points  should  be  chosen  to  represent  adequately  the  behaviour  of  the  curve.  If 
the  data  points  are  provided  they  must  be  supplemented  from  a hand  drawn  curve 
if  there  are  appreciably  less  than  20  of  them.  The  shape  of  the  curve  gives 
some  guidance  on  the  placing  of  the  interior  knots.  There  should  not  be  too 
many  of  these  as  they  can  produce  unwanted  fluctuations.  The  number  of  data 
points  must  exceed  the  number  of  interior  knots  by  at  least  A.  In  general,  more 
knots  are  needed  where  the  behaviour  of  the  curve  is  changing  rapidly,  and  less 
where  it  is  changing  slowly.  The  exact  positioning  of  the  knots  is  not  critical. 
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As  the  program  ODNOZ  calculates  the  axial  position  x from  the  area  A 
(when  performing  the  throat  calculation)  and  finds  A from  x at  stations 
downstream  of  the  throat,  two  cubic  spline  functions  r * s^(x)  and  x - sr(r) 
must  be  fitted  to  the  set  of  data  points  and  hence  two  sets  of  knots,  at  given 
axial  locations  for  the  function  in  x , and  at  given  radial  points  for  the 
function  in  r , must  be  provided. 

An  example  of  a contoured  nozzle,  providing  guidance  on  the  placing  of 
interior  knots,  is  given  in  Section  6. 

ODNOZ  and  TRANOZ  each  occupy  about  30  K words  of  store  (decimal)  and 
SUPERNOZ  requires  1*0  K.  Typical  run  times  for  the  programs  on  an  I CL  1904S 
computer  are  one  minute  for  ODNOZ  and  two  minutes  for  TRANOZ.  The  time  taken  to 
run  SUPERNOZ  with  no  chemical  reactions  (frozen  chemistry)  is  about  one  minute, 
but  when  chemical  reactions  are  included  the  time  varies  considerably,  depending 
on  the  number  of  species  and  reactions.  A typical  run  time  is  probably  about 
one  hour. 

5. I ODNOZ 

ODNOZ  has  two  input  channels  (5  and  7)  and  two  output  channels  (6  and  9). 
Channel  5 is  used  to  input  the  basic  data  and  channel  7 is  used  to  input 
THERMDATA.  Channel  6 produces  the  general  ( 1 i nep r i n te r)  output,  comprising 
chamber  and  throat  conditions,  the  solution  at  a specified  number  of  axial 
stations  downstream  of  the  throat  and  the  equilibrium  exit  plane  conditions. 
Channel  9 produces  the  table  of  flow  properties,  starting  at  the  throat  (M»l), 
required  for  input  to  TRANOZ.  Care  must  be  taken  to  ensure  that  the  range  of 
Much  numbers  in  the  table  spans  the  value  likely  to  be  encountered  at  the  nozzle 
throat  wall  point.  This  could  range  from  slightly  greater  than  unity  for  a 
nozzle  with  a large  normalised  throat  wall  radius  of  curvature  (R*/r*  > 3,  say), 
to  about  1.3  for  a nozzle  with  a small  normalised  throat  wall  radius  of 
curvature  (R'/r*  •»  0.5,  say). 

INPUT  CHANNEL  9 

Input 

No. 


1 

FORMAT 

IHEAD 

(12A6) 

run  title. 

2 

FORMAT 

(510) 

1 AX 

2 

number  of  axial  stations  at  which  the  solution 

is  required  to  construct  the  table  of  flow 
properties  for  input  to  TRANOZ  (max.  20). 
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1 nput 
No. 

ICONT 

0 for  conical  nozzle,  1 for  contoured  nozzle. 

ISTOT 

J 

total  number  of  chemical  species  (max.  ^0). 

1C 

: 

number  of  fundamental  species  excluding  those 
to  be  considered  chemically  inert. 

INERT 

: 

number  of  inert  species. 

3 

FORMAT 

(12E0.0) 

TOD 

initial  guess  at  chamber  temperature,  K. 

POD 

chamber  pressure,  atmospheres. 

HIN 

specific  enthalpy  of  the  propellant,  cal  g \ 

RL 

throat  radius,  cm. 

XSTART 

starting  point,  downstream  of  the  throat,  at 
which  the  solution  is  required  for  construc- 
tion of  the  table  of  flow  properties  (usually 
0.02  r*) , cm. 

XFIN 

'• 

distance  from  throat  to  nozzle  exit  plane,  cm. 

XINC 

x increment  between  points  in  the  table  of 
flow  properties  (usually  0.02  r*) , cm. 

ACC 

: 

accuracy,  usually  1.0  E-06. 

NOZZLE 

GEOMETRY 

(see  Figs 

5a  and  5b) 

FORMAT 

(12E0.0) 

TH 1 

first  inlet  cone  half  angle,  ueg. 

XXI 

•* 

distance  from  throat  (-ve)  at  which  first  and 
second  inlet  cones  meet,  cm. 

TH2 

•’ 

second  inlet  cone  half  angle,  deg. 

X2 

distance  from  throat  (-ve)  at  which  second 
inlet  cone  meets  first  circular  arc  throat 
section,  cm. 

RADI 

♦ 

first  circular  arc  radius,  cm. 

X3 

: 

distance  from  throat  at  which  first  and  second 
circular  arcs  meet,  cm.  (+ve,  -ve  or  zero). 

Input 

is  incl uded  only 

if  the  nozzle  has  a conical  exit  (ICONT  = 0) . 
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Input 

No. 

5 

FORMAT  (12E0.0) 

RAD2  : 

second  circular  arc  radius,  cm. 

X*4  ; 

distance  from  throat  to  end  of  second  circular 
arc  section,  cm. 

TH3  : 

first  exit  cone  half  angle,  deg. 

' 

X5  : 

distance  from  throat  at  which  the  first  and 
second  exit  cones  meet,  cm. 

j 

TH*4  : 

second  exit  cone  half  angle,  deg. 

1 npu  t 

group  6 is  included  only  i 

f the  nozzle  has  a contoured  exit  (ICONT  = l). 

6.1 

FORMAT  (510) 

i 

NP  : 

number  of  data  points  (max.  *40). 

4 

i 

t 

NKX  : 

number  of  intervals  ( = number  of  interior 
knots  + 1)  for  the  cubic  spline  approximation 
r = s^(x)  (max.  8) . 

NKR  : 

number  of  intervals  for  the  cubic  spline 
approximation  x = sr(r)  (max.  8). 

6.2 

FORMAT  (12E0.0) 

! 

There  now  follow  NP  records 

containing  the  set  of  data  points  (x.  , r.) 

wi  th 

weights  w. 

1 

XC  ( 1 ) : 

axial  co-ordinate  of  the  i th  data  point,  cm. 

Q(D  : 

radial  co-ordinate  of  the  i th  data  point,  cm. 

WT(l)  : 

weight  of  the  i th  data  point  (usually 

Wj=1.0  , but  if  the  curve  is  required  to 

pass  exactly  through  a particular  point,  e.g. 

throat  or  exit  plane,  set  w.  > 1.0  ). 

for  1=1 , 2 , ...  , NP 

6.3  FORMAT  (12EO.O) 


The  next  record  contains  the  NKX-1  internal  knots  for  fitting  the  spline 

function  r = s (x) 
x 


RKX(5) 

RKX(J) 

RKX (NKX+3) 


axial  position  of  the  first  internal  knot,  cm. 
axial  position  of  the  i th  internal  knot,  cm. 
axial  position  at  the  last  internal  knot,  cm. 


2*4 
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I nput 
No. 

Eight  additional  knots  are  inserted  by  the  program,  the  first  four  are  set 

equal  to  *,  , the  last  four  to  x 

1 np 

6.4  FORMAT  (12E0.0) 

The  next  record  contains  the  NKR-1  internal  knots  for  fitting  the  spline 
func  tion  x * sr( r) 

RKR(5)  : radial  position  of  the  first  internal  knot,  cm. 

RKR(J)  : radial  position  of  the  jth  internal  knot,  cm. 

RKR(NKR+3) : radial  position  of  the  last  internal  knot,  cm. 

Eight  additional  knots  are  inserted  by  the  program;  the  first  four  are  set 

equal  to  r,  , the  last  four  to  r 

1 np 

Input  group  7 contains  the  chemical  data. 

7.1  FORMAT  (A6.E12.6) 


There  now  follow  1C  records  containing  the  names  and  concentrations  of  the 
fundamental  species,  and  INERT  records  containing  the  names  and  concentrations 
of  the  inert  species. 


I AT  I T ( I ) 


name  of  the  i th  fundamental  species. 


Q(  I ) 


number  of  gram  atoms  of  the  i th  fundamental 
species/100  grams  of  propellant. 


1 = 1,2 1C 

I AT  I T ( I ) : 

Q(  I ) : 


name  of  the  i th  (inert)  species. 

number  of  gram  atoms  of  the  i th  (inert) 
species/100  grams  of  propellant. 


1= I C + 1 , ...  , 1C  + INERT 
7.2  F0RMAT(A6 , 1 8F4 . 1 ) 


The  remaining  I STOT- I C- INERT  records  contain  the  names  of  the  dependent 
species,  along  with  the  matrix  A , which  describes  the  formation  of  the  i th 
species  in  terms  of  the  jth  fundamental  species. 
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I AT  I T ( I ) 


name  of  the  i th  (dependent)  species. 


ALPHA(I ,1) : 


number  of  moles  of  the  first  independent 
species  per  mole  of  the  i th  species. 


ALPHA( I ,J) : 


number  of  moles  of  the  j th  independent  species 
per  mole  of  the  i th  species. 


ALPHA ( I , 1C) : 

for  I = I C+ 1 NE  RT+ 1 , ISTOT 


number  of  moles  of  the  last  independent 
species  per  mole  of  the  i th  species. 


INPUT  CHANNEL  7 

Thermodynamic  data  for  chemical  species  likely  to  be  required  in  a rocket 
nozzle  calculation  are  stored  in  a standard  data  file,  in  blocks  of  A records  in 
the  form  specified  in  Ref.  10.  The  thermodynamic  data  for  the  species  in  the 
calculation  are  selected  by  comparing  the  6 character  name  of  the  first  block  in 
the  standard  data  file  with  each  of  the  species  names  previously  specified.  If 
there  is  no  match  the  block  of  data  is  discarded  and  the  process  is  repeated 
with  the  next  block  until  the  thermodynamic  data  for  all  the  species  present  in 
the  calculation  have  been  found.  If  the  calculation  includes  a species  not 
appearing  in  the  standard  thermodynamic  data  file,  or  if  the  name  is  different, 
the  message 

"(species  name)  not  in  thermodynamic  data  file" 
is  output  and  the  calculation  is  terminated. 

The  required  thermodynamic  data  are  input  in  the  form: 

Input 

No, 

1 F0RMAT(A6) 

IN  : name  of  species  in  standard  thermodynamic 

data  file. 

The  following  records  are  saved  only  if  I N= I AT  I T ( J ) 


2 F0RMAT(4E 16.8) 
WT  ( J ) ; 

AC (9 , J)  : 


molecular  weight  of  jth  species, 
coefficient  a^  in  polynomial  expression 
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I nput 
No. 


CS(J)  : 

CH  ( J)  : 

I 

3 F0RMAT(4E 16 . 8) 

AC ( I , J)  : 

1*=)  ,4 

l*  FORMAT  ( ^ t 16.8) 

AC(I,J)  : 

1-5,8 


5.2  TRANQZ 

TRANOZ  has  two  input  channels  (1  and  3)  and  two  output  channels  (2  and  4). 
Channel  I is  used  to  input  the  data  defining  the  transonic  flow  field  and  the 
axial  location  of  the  point  on  the  nozzle  wall  from  which  the  initial  data  line 
is  constructed  for  starting  the  method  of  characteristics  solution  (see  Fig.  2). 
The  point  chosen  to  start  the  supersonic  solution  is  usually  the  throat  wall 
point  since  the  Mach  number  here  is  always  greater  than  unity  for  particle-free 
flow,  and  this  choice  permits  a finite  rote  chemical  treatment  as  far  upstream 
as  possible.  The  transonic  flow  domain  must  be  large  enough  to  contain  the 
whole  of  the  constant  Mach  number  method  of  characteristics  starting  line,  but 
must  not  extend  too  far  downstream  because  for  large  Mach  numbers  the  solution 
procedure  becomes  highly  unstable  and  requires  a very  small  relaxation  factor, 
and  hence  many  iterations,  to  obtain  convergence.  If  the  upstream  and  down- 
stream boundaries  are  placed  at  x * -0.6  r*  and  x - 0.8  r,v  respectively, 
both  these  requirements  should  be  satisfied  for  all  nozzles  likely  to  be 
encountered,  with  gas  phase  flow.  A relaxation  factor  of  0.3  is  optimum  for 
this  configuration,  convergence  being  achieved  in  all  cases  in  far  fewer  than 
100  iterations.  The  maximum  number  of  iterations  to  be  permitted  in  the 


4 

CP(J)  - 

for  specific  heat  of  jth  species  (cal  mol 

k'V 

thermodynamic  property  (cal  mol  ' K 

thermodynamic  property  C (cal  nx)l  ' K 

coefficients  a_^  to  a_ ^ in  polynomial 
expression  for  specific  heat  of  jth  species 
(cal  mol  ' K n) . 

coefficient  a^  to  a^  in  polynomial 
expression  for  specific  heat  of  jth  species 
(cal  mol  ' K n) . 
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calculation  and  the  interval  between  iterations  at  which  the  solution  is  to  be 
printed  are  also  input. 

Channel  3 is  used  to  input  the  data  prepared  by  ODNOZ.  Channel  2 outputs 
the  whole  transonic  flow  field  solution,  giving  either  the  final  solution  only, 
or  including  the  solution  at  intermediate  iterations  if  required.  If  the 
iterative  process  does  not  converge  to  the  required  accuracy,  the  message 

“maximum  number  of  iterations  reached" 

is  output  and  the  solution  at  that  stage  is  given.  The  calculation  may  also  be 
terminated  with  the  message 

“no  solution  to  mass  flux  equation  at  x = “ 

which  means  that  the  calculated  mass  flux  at  a station  just  downstream  of  the 
throat  is  less  than  the  required  mass  flux.  This  may  happen  with  nozzles  which 

. 

have  a very  small  normalised  throat  wall  radius  of  curvature  (R-'/r*  < l)  , since 
in  the  first  iteration  the  streamline  minima  are  assumed  to  be  at  the  throat, 
and  for  such  nozzles  they  are  often  well  downstream  of  the  throat.  Hence,  at  an 
axial  station  just  downstream  of  the  throat,  some  of  the  streamlines  may  have 
negative  gradients  rather  than  positive  ones,  thus  causing  a large  error  in  the 
mass  flux  integration.  For  this  reason  the  axial  step  size  if  always  increased 
to  0.2  r*  if  the  user  sets  it  to  less  than  that.  The  program  may  be  restarted 
by  increasing  the  nozzle  step  size  to  slightly  greater  than  0.2  r*,  to  by-pass 
the  region  where  the  streamlines  are  still  converging. 

Channel  **  produces  the  data  required  for  input  to  SUPERNOZ,  comprising 
chamber  conditions,  thermochemical  data  and  the  properties  along  the  constant 
Mach  number  line.  A check  must  be  made  to  ensure  that  the  Mach  number  at  the 
point  on  the  nozzle  wall  chosen  to  start  the  supersonic  solution  is  contained 
within  the  table  produced  by  ODNOZ.  Care  must  also  be  taken  to  ensure  that  the 
whole  of  the  constant  Mach  number  line  lies  within  the  calculated  transonic 
domain.  For  nozzles  with  a small  normalised  throat  wall  radius  of  curvature  the 
point  where  the  line  crosses  the  axis  can  be  well  downstream  of  the  point  where 
it  crosses  the  wall  (see,  for  example.  Fig.  2). 
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INPUT  CHANNEL  1 

I nput 
No. 


FORMAT (I00E0.0) 

XSTART  : point  on  nozzle  wall  from  which  the  initial 

data  line  is  constructed  for  starting  the 
supersonic  solution  (usually  0.0),  cm. 


XSTEP  : x increment  between  axial  stations  in 


transonic  flow  field  (usually  0.2  r* , If 
xstep  < 0.2  r*  , program  sets  xstep  - 
0.2  r*  ),  cm. 

XM1N  : 

minimum  value  of  x in  transonic  flow  field 
(usually  -0.6  r*) , cm. 

XMAX  : 

maximum  value  of  x in  transonic  flow  field 
(usual  1 y 0.8  rA) , cm. 

RF  : 

relaxation  factor  (usually  0.3). 

ACC  : 

accuracy,  used  to  test  the  convergence  of  the 
Mach  number,  usually  1.0  E-05. 

FORMAT (310) 

NN  : 

number  of  radial  points  in  the  transonic  and 
supersonic  flow  fields  (maximum  16) . 

ITER  : 

maximum  number  of  iterations  to  be  permitted 
in  the  calculation. 

1 PR  1 NT  : 

number  of  iterations  between  print  stations. 

5.3  SUPERNOZ 

SUPERNOZ  has  three  input  channels  (1,  5 and  7)  and  two  output  channels 
(2  and  4) . Channel  1 is  used  to  input  the  data  produced  by  TRANOZ,  channel  5 is 
used  to  input  the  number  of  chemical  reactions  and  the  name  of  each  reaction  as 
it  appears  in  the  standard  reaction  data  file  REACTIONS,  and  channel  7 is  used 
to  input  the  standard  reaction  data  file,  containing  the  reaction  mechanisms  and 
rate  coefficients  of  most  reactions  likely  to  be  significant  inside  rocket 
nozzles''.  The  standard  file  is  searched  by  the  program  to  find  those  reactions 
to  be  included  in  the  calculation. 


Channel  2 produces  the  general  output,  comprising  the  non-equilibrium 
solution  along  each  data  line  in  the  flow  (see  Fig.  4)  and  the  non-equilibrium, 
two-dimensional  exit  plane  conditions.  Channel  4 outputs  the  exit  plane 
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conditions  for  use  in  preparing  the  input  to  a rocket  exhaust  plume  structure 

I 2 

program,  such  as  REP  or  8AFL  . 

INPUT  CHANNEL  5 

Input 

No. 

1 FORMAT (10) 

NR  : total  number  of  chemical  reactions  to  be 

included  in  the  calculation  (max  40). 

If  NR-0  the  gas  composition  is  frozen,  and  no  reaction  data  are  required. 

2 FORMAT ( I 3A6) 

NUMR(l)  : number  of  reaction  (e.g.  A1 , A2 , ...) 

as  It  appears  In  the  standard  reaction  data 
1-1 , NR  f| le  and  in  Ref.  1 1 . 

INPUT  CHANNEL  7 

The  standard  reaction  data  file  REACTIONS  contains  all  the  reactions,  with 
rate  coefficients,  which  are  listed  In  Ref.  II.  An  example  of  a reaction 
rivchanism  with  rate  coefficients  is  yfven  in  Table  2.  The  reaction  data  are 
input  as 


FORMAT  (A6 ,2  X,6  (Afe , 1 X) , IX,  E8.2,  F4.1,  F9.1,  3X.  12) 


NUM 

NAME (J) 
J-1,3 


number  of  reaction  (e.g.  A1,A2,  ...  etc.) 

names  of  first,  second  and  third 
species  on  the  left  side  of  the  reaction. 


1 1 


NAME (j) 
J-4,6 


names  of  first,  second  and  third 
species  on  the  right  side  of  the  reaction. 


AC1 


AC  2 

AC  3 


pre-exponential  function  R in  the  reaction 
rate  expression  k^.  - 0.  Tn  exp(E/T) 

(mol ecul e-cm-g  units). 

the  exponent  a in  the  above  reaction  rate 
expression  (dimensionless). 

the  Arrhenius  factor  E In  the  above 
expression  (K) . 


N 1 


1 for  two-body  reaction,  2 for  three-body 
reaction. 
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If  NUM  is  not  one  of  the  reactions  to  be  included  in  the  calculation,  the 
data  are  discarded  and  the  next  entry  in  the  standard  reaction  data  file  is 
input.  However,  if  NUM  » NUMR(l),  the  names  of  the  species  on  the  left  and 
right  hand  sides  of  the  reaction  equation  are  compared  with  the  names  of  the 

and 


species  in  the  calculation  and  the  stoichiometric  coefficients  v.  . 

in  the  reaction  i on  the  reactants  and  products  sides 


v:  . 
' J 


of  species  j 


respectively,  are 

set 

to 

1 or  zero  accord 

ingly.  The  reaction  rate  coefficient 

data  are  saved  as 

» 

j 

K1 

(1) 

«*  AC  1 

23 

if  N1  = 1 or 

i 

K 1 

(1) 

= AC  1 

x 6|0  23 

x 6]0  23  if  N1  = 2 , 

K2 

(1) 

= AC2 

and  K3(l) 

= AC3  . 

Then 


K1  TK^  exp(«3/T)  cc  mol  ^ sec  ' 


or 


2 .-2  -1 
cc  mol  sec 


If  any  of  the  reactions  to  be  included  in  the  calculation  cannot  be  found  in  the 
reaction  data  file,  the  message 

"(reaction)  not  found  in  reaction  file" 

is  output  and  the  calculation  is  terminated. 

6 SAMPLE  CALCULATION 


A sample  calculation  was  carried  out  to  find  the  two-dimensional,  non- 
equilibrium exit  plane  conditions  of  a liquid  bipropellant  engine.  Table  1 
contains  the  propellant  composition  and  chamber  pressure.  Table  2 lists  the 
chemical  reactions  used  in  the  supersonic  calculation,  and  the  nozzle  geometry 
is  shown  in  Fig.  6a.  The  nozzle  has  a contoured  exit  section,  along  which  a 
table  of  (x  , r)  values  are  given.  The  data  points,  and  the  selected  axial  and 
radial  knots,  are  shown  in  Fig.  6b.  The  input  data  for  the  program  CDNOZ,  which 
include  the  names  and  concentrations  of  the  fundamental  species  and  the  names 
and  composition  (in  terms  of  the  fundamental  species)  of  the  dependent  species 
to  be  included  in  the  calculation  are  shown  in  Fig.  7.  A typical  output  station, 
showing  the  one-dimensional  equilibrium  throat  conditions,  is  given  in  Fig.  8. 

The  transonic  flow  field,  with  contours  of  constant  Mach  number,  is  shown  in 
Fig.  9.  The  nozzle  of  the  test  case  engine  has  a normalised  throat  wall  radius 
of  curvature  of  1.38,  which  is  relatively  small,  but  the  constant  Mach  number 
line  emanating  from  the  throat  wall  point  (M  - 1.17),  chosen  to  start  the  method 
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of  characteristics  solution,  is  shown  in  Fig.  9 to  lie  well  within  the  transonic 
solution  domain. 

The  two-di  nvns ional  exit  plane  profiles  of  pressure  and  temperature  produced 

by  SUPERNOZ  are  shown  in  Figs  10  and  11  respectively,  compared  with  the  one- 

dimensional  results  from  FIRAC^.  The  value  of  a two-dimensional  treatment  is 

immediately  apparent.  The  one-dimensional  program  calculates  the  exit  plane 

-2 

pressure  to  be  constant  at  0.86  5 Nm  , whereas  SUPERNOZ  calculates  it  to  lie 

-2  U -2 
between  0.45^  5 Nm  on  the  nozzle  axis  and  2.151Q  5 Nm  on  the  wall. 

Similarly,  the  one-d i mens iona 1 exit  plane  temperature  is  constant  at  1790  K 

whilst  that  from  the  two-dimensional  program  varies  between  1680  K and  2250  K. 

The  variations  of  the  gas  properties  along  the  nozzle  axis  and  wall  are 
illustrated  in  Figs  12  to  14.  Fig.  12  shows  the  axis  and  wall  distribution  of 
pressure.  The  axis  pressure  drops  at  a fairly  constant  rate  throughout  the 
nozzle,  whereas  the  wall  pressure  falls  initially  very  rapidly  in  the  convergent 
part  of  the  nozzle  and  then  much  more  slowly  in  the  expansion  region  until,  when 
the  compression  corner  is  reached,  it  begins  to  increase.  Similar  fluctuation 
is  exhibited  by  the  temperature,  as  shown  in  Fig.  13.  The  variation  of  the  gas 
velocity  along  the  nozzle  is  shown  in  Fig.  I1*,  from  which  it  can  also  be  seen 
that  the  exit  velocity  varies  between  2845  ms  1 on  the  axis  and  2415  ms  1 on 
the  wall. 

The  engine  thrust  calculated  by  the  two-dimensional  program  is  not  very 
different  from  that  calculated  with  the  one-dimensional  program.  SUPERNOZ 
computed  the  thrust  to  be  2110  N,  whilst  FIRAC^  gave  it  as  2090  N.  However,  the 
exhaust  plume  of  the  test  case  engine  calculated  by  the  REP  program1  with  two- 
dimensional  nozzle  exit  plane  conditions  as  input  differs  significantly  from  that 
calculated  with  uniform  exit  plane  conditions  as  input.  Fig.  15  shows  a 
comparison  of  the  predicted  plume  centreline  temperature  profile  for  the  engine 
under  sea  level,  static  conditions  with  both  uniform  and  non-uniform  nozzle  exit 
plane  conditions  as  input.  Downstreamof  the  shock  region  the  temperature 
profiles  are  similar,  but  close  to  the  nozzle  exit  plane  they  are  very  different. 
Although  the  number  and  locations  of  shock  waves  calculated  by  the  program  with 
the  two  different  sets  of  inlet  conditions  are  the  same,  the  strengths  of  the 
shocks  calculated  with  the  two-dimensional  exit  plane  conditions  as  input  are 
much  greater.  This  could  be  quite  significant  in  a number  of  contexts,  particu- 
larly in  a fusing  context,  where  accurate  knowledge  of  the  strengths  and 
locations  of  temperature  discontinuities  is  very  important. 
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7 CONCLUSIONS 

A set  of  computer  programs  has  been  written  to  calculate  the  two-dimensional 
flow  conditions  in  axisymmetric  nozzles.  The  chemical  species  are  assumed  to  be 
in  equilibrium  in  the  combustion  chamber  and  as  far  as  the  nozzle  throat,  whilst 
frozen  or  finite  rate  chemistry  may  be  applied  in  the  supersonic  region  of  the 
nozzle.  A summary  of  the  theoretical  basis  of  the  programs  is  included,  with 
full  operating  instructions.  A sample  calculation  is  presented,  illustrating 
the  wide  radial  variation  in  rocket  exhaust  properties  which  may  occur  at  the 
nozzle  exi t plane. 


Although  the  programs  do  not  produce  a significant  improvement  over  previous 

methods^  in  accuracy  of  thrust  calculations,  they  provide  far  more  realistic 

estimates  of  nozzle  exit  plane  properties  and  hence  better  starting  conditions 

1 2 

for  an  exhaust  plume  structure  analysis  program  ' . They  are  particularly 
important  in  helping  to  provide  accurate  predictions  of  the  shock  structure  of 
a rocket  exhaust  plume. 
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TABLE  1 

Test  case  parameters 

Fue  I 

Mono-methyl  hydrazine  ^H6N2 

Oxi dant 

Dinitrogen  tetroxide 

1.8 

3.^5  MPa  (3**  atm) 
1.139  cm 
1.575  cm 
1.38 

2.808  cm 


Oxi di ser/fuel  ratio 

Chamber  pressure 

Nozzle  throat  radius 

Nozzle  throat  radius  of  curvature 

R*/r* 

Nozzle  exi t rad i us 


35 

UNLIMITED 


UNLIMITED 


' 


1 

i 

i 

i 

! 

: 

I 

< 

i 


: 


i 


TABLE  2 

Reaction  mechanism  for  test  case 


M 11 

No . 

Reaction 

Rate  coefficientu 

Uncertainty 

factor 

A1 

0 + 0 + M 

-> 

0,  -*•  M 

3 x 

10"3A  exp(900/T) 

10 

A 2 

0 + H + M 

-► 

OH  + M 

1 X 

10-29  T-1 

30 

AT 

H + H ♦ M 

-► 

H„  + M 

3 x 

10“30  T'1 

30 

A4 

H + OH  ♦ M 

-► 

z 

H00  -*■  M 

2 x 

-25  -2 

10  T 

10 

AS 

CO  ♦ 0 ♦ M 

z 

CO.,  + M 

7 x 

10”33  exp(-2200/T) 

30 

Bl 

OH  H2 

-> 

H20  + H 

3.6 

x 10_U  exp(-2600/T) 

2 

B2 

0 + H., 
H + 02 

-> 

OH  + H 

3.0 

x io'14  T exp(-4480/T) 

1.5 

B 1 

-> 

OH  + 0 

3.7 

x 10~10  exp(-8450/T) 

1.5 

B4 

CO  + OH 

-> 

CO.,  + H 

2.8 

x 10“17  T1*3  exp(330/T) 

3 

B5 

OH  + OH 

h2o  + 0 

1 x 

10_U  exp(-550/T) 

3 

PI 

H ♦ 0.,  + M 

> 

H02  + M 

2 x 

10"32  exp(500/T) 

10 

P3 

H + H02 

-> 

OH  ♦ OH 

4 x 

10“10  exp(-950/T) 

5 

P4 

H + H02 

-► 

»2  + °2 

4 x 

10-U  exp(-350/T) 

5 

P5 

h2  ♦ ho2 

-► 

»20  ♦ OH 

1 x 

10-12  exp(-9400/T) 

10 

P6 

co  + ho2 
o ho2 

- V 

CO  + OH 
z 

2.5 

x 10_1°  exp(-11900/T) 

10 

P7 

-> 

OH  + O0 

8 x 

10-U  exp(-500/T) 

30 

P8 

OH  ♦ HO., 

-► 

h2o  + 02 

1 x 

10_1° 

30 

Q1 

CHO  ♦ 0 

-> 

CO  + OH 

3 x 

10-11  exp(-2500/T) 

100 

Q2 

CHO  + OH 

-► 

CO  ♦ H.,0 

1 x 

io-11 

100 

Q3 

CO  ♦ H ♦ M 

z 

CHO  + M 

2 x 

10~33  exp(-850/T) 

100 

Q4 

CHO  ♦ H 

-► 

CO  + h2 

2 x 

10-10  exp(-2500/T) 

30 

U3 

NO  + H ♦ M 

-► 

HNO  + M 

1.5 

x 10~32  exp(300/T) 

10 

U4 

HNO  + H 

-► 

NO  ♦ H, 

8 x 

io"12 

10 

U5 

HNO  + OH 

-► 

z 

NO  + H.,0 

6 x 

io-11 

5 

U6 

N-  ♦ 0 

-> 

z 

NO  + N 

1.3 

x 10~10  exp(-38000/T) 

5 

U7 

N ♦ 02 

-> 

NO  ♦ 0 

1.1 

x 10~14  T exp(-3150/T) 

3 

* 


a ml-molecule-second  units 

b defined  as  in  Ref.  II 
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uni  iMiun 


Numeric  I aturc 


nozzle  cross-sectional  area  normalised  by  the  throat  area;  diabatic  heat 
addition  term  linking  fluid  dynamic  and  relaxation  processes 

energy  exchange  term  linking  fluid  dynamic  and  relaxation  processes 

number  of  independent  chemical  species;  speed  of  sound,  ms  1 

mass  fraction  of  jth  species 

speci f i c heat , J kg  1 


derivative  with  respect  to  x of  Y.  , dis- 
function defined  by  equation  (14.17) 
function  defined  by  equation  ( 4 . 1 8) 
enthalpy  of  the  jth  species,  J mol  1 

total  enthalpy,  J mol  1 ; function  defined  by  equation  (14.19) 

forward  rate  coefficient  for  the  i th  reaction,  cc  mol  1 s 1 or 

2 .-2  -1 
cc  mol  s 

equilibrium  constant  for  the  jth  species,  based  on  partial  pressures 
equilibrium  constant  of  the  i th  reaction 
mixture  molecular  weight,  g mol  1 
molecular  weight  of  jth  species,  g mol  ^ 

Mach  number 

initial  mole  fraction  of  species  j 

initial  number  of  moles  of  jth  species  per  unit  mass  of  propellant 
pressure,  Nm  ^ 
pressure,  atmospheres 

number  of  chemical  reactions;  radial  co-ordinate,  normalised  by  throat 
radius 

throat  radius,  cm 

universal  gas  constant,  8.3143  J K ' mol  ' , 82.  056  cc  atm  mol  1 K 1 
nozzle  throat  wall  radius  of  curvature 


total  number  of  chemical  species 
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UNLIMITED 


UNLIMITED 


entropy  ot  jth  species  at  1 atmosphere,  cal  mol  * K ' 
mixture  entropy,  non-dimensional 
temperature,  K 

non-dimensional  temperature,  ° T/Tq 

. . , . -1 
axial  velocity,  ms 

radial  velocity,  ms  1 

, • "I 

velocity,  ms 

net  species  production  rate 

axial  distance  from  nozzle  throat,  normalised  by  throat  radius 
mole  fraction  of  species  j 

Mach  angle,  angle  between  streamline  and  Mach  line  characteristics 
number  of  moles  of  species  j per  mole  of  species  i 

difference  between  the  stoichiometric  coefficients  of  the  jth  species 

on  the  product  and  reactant  sides  of  the  i th  reaction,  B. . = v!  . - v.  . 

|J  'J  'J 

j-i 

isentropic  exponent 

concentration  of  jth  species,  mol  g ' 

Kronecher  delta 
dens i ty , kg  m ^ 
s t reaml  i ne  angl  e 

chemical  potential  of  jth  species  evaluated  under  standard  pressure 
conditions  (1  atmosphere) 

stoichiometric  coefficient  of  the  jth  species  on  the  reactant  side  of 
reaction  i 

stoichiometric  coefficient  of  the  jth  species  on  the  product  side  of 
reaction  i 


E«.i 
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nozzle  wall. 


COMBUSTION 


EXHAUST 


CHAMBER 


NOZZLE 


PLUMt 


throat 


solution  field  conditions  properties 

FIG.1  CONVERGENT/DIVERGENT  NOZZLE  WITH  CALCULATION  SEQUENCE 


nozzle  axis 


FIG.  2 LOCATIONS  OF  SONIC  LINE  AND  METHOD  OF  CHARACTERISTICS 
STARTING  LINE 


FIG.3&4 


FIG.  3 FLOW  FIELD  CONFIGURATION  FOR  TRANSONIC  SOLUTION 


FIG.  *4  CHARACTERISTIC  MESH  FOR  SUPERSONIC  SOLUTION 


Radius,  cm 


FIG.6A&B 


Oistance  from  nozzle  throaf,  cm 


FIG.6A  NOZZLE  GEOMETRY  OF  TEST  CASE  ENGINE 


FIG.6B  NOZZLE  CONTOUR  SHOWING  PLACEMENT  OF  KNOTS 
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FIG.  7 


w 


TEST  CASE 
20  1 14  4 0 

3200.0  34.0229953  67.8144  1.  13919  0.  08  5. 
33.233  -1.309144  33.233  -1.309144  1.5748 
12  4 5 


50  672 

0.0 


0.02  1 


0.0  1.13919 

0.42672  1.2065 
0.93472  1.41732 
1.44272  1.65227 
1.95072  1.88468 
2.45872  2.09677 


5.0 

1.0 

1 l 

I i 
1 » 

1 I 


0 
0 
0 
0 

2.96672  2.28727  1.0 
3.47472  2.45237  1.0 
3.98272  2.58953  1.0 
4.49072  2.69875  1.0 
4.99872  2.77495  1.0 
5.50672  2.80797  5.0 
0.93472  2.45872  4.49072 
1.4  1 732  1.8  8 4 68  2.4  5 2 3 7 


> CONTOUR  DATA 


2.  69875 


CO 

7. 7509 1 2E-0 1 

H2 

2.325273 

02 

1*009744 

N2 

1.473736 

C02 

1.0  0.0 

0.5 

0.0 

OH 

0.0  0.5 

0.5 

0.0 

H 

0.0  0.5 

0.0 

0.0 

O 

0.0  0.0 

0.5 

0.0 

H02 

0.0  0.5 

1.0 

0.0 

CHO 

1.0  0.5 

0.0 

0.0 

H20 

0.0  1.0 

0.5 

0.0 

N 

0.0  0.0 

0.0 

0.5 

NO 

0.0  0.0 

0.5 

0.5 

HNO 

0.0  0.5 

0.5 

0.5 

FUNDAMENTAL  SPECIES 


> DEPENDENT  SPECIES 


FIG.  7 


INPUT  DATA  FOR  PROGRAM  ODNOZ 


. 0E-6 


» 


EQUILIBRIUM  THROAT  CONDITIONS 


TITLE*-  TEST  CASE 
ACCURACY  1.0QE-Q6 


MOLE  FRACTION  NO. PER  C.C. 


CO 

1 . 182430E-0 1 

5.  550612E 

18 

H2 

1.  199625E-0  1 

5. 631327E 

18 

02 

1. 1I6751E-03 

5.  242298E 

16 

N2 

3.  1 56Q99E-0  1 

1 • 48 1 549 E 

19 

C02 

4.833789E-02 

2.  2 6909 7E 

18 

OH 

1.4110 50 E-02 

6. 623807E 

17 

H 

1.321 660  E-Q  2 

6.204190E 

17 

0 

9. 18299 1 E-0  4 

4.310717E 

16 

H02 

1.852384E-06 

8.  695536E 

13 

CHO 

1. 68  9 2 4 0 E-0  5 

7.929  700E 

14 

H20 

3. 661 572E-0 1 

1 . 718830E 

19 

N 

1 .8  29355E-0  6 

8. 587432E 

13 

NO 

2.30 543 6E-0  3 

l. 082227E 

17 

HNO 

1.471 348 E-0  6 

6.9  0 68  62E 

13 

TEMPERATURE 

3029.93 

K 

ENTHALPY 

-8451.1 

J 

MQL-1 

MOLECULAR  VT. 

21.4789 

G MOL- 1 

ENTROPY 

11.99101 

J 

MOL- 1 K-l 

SPECIFIC  HEAT 

2099.33 

J KG-  1 

DENSITY 

1.67394 

KG  M-3 

VELOCITY 

1 1 67. 01 

M 

SEC-1 

PRESSURE 

19.3763 

ATM 

MACH  NO. 

1.00730 

GAMMA 

1.13777 

CSTAR 

1741.63 

M SEC-1 

MASS  FLOW 

0.796 

KG  SEC-1 

FIG.  8 TYPICAL  OUTPUT  STATION  FROM  PROGRAM  ODNOZ  SHOWING  ONE- 
DIMENSIONAL EQUILIBRIUM  THROAT  CONDITIONS 
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Abstract 

This  report  describes  a set  of  computer  programs  which  calculates  the  inviscid 
two-dimensional  expansion  of  gaseous  exhaust  mixtures  through  rocket  nozzles. 
Equilibrium  chemistry  in  the  combustion  chamber  and  as  far  as  the  nozzle  throat 
is  assumed,  whilst  finite  rate  or  frozen  chemistry  may  be  applied  from  the  throat 
onwards.  The  transonic  flow  field  is  calculated  by  the  method  of  streamline 
curvature  and  the  flow  In  the  supersonic  region  of  the  nozzle  is  solved  by  the 
method  of  characteristics. 


The  theory  on  which  the  programs  are  based  is  summarised  and  a comprehensive 
user's  guide  is  included.  A sample  calculation  is  presented. 


